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Abstract. We model the process of dust coagulation in protoplanetary disks and calculate how it affects their 
observational appearance. Our model involves the detailed solution of the coagulation equation at every location 
in the disk. At regular time intervals we feed the resulting 3-D dust distribution functions into a continuum 
radiative transfer code to obtain spectral energy distributions. We find that, even if only the very basic - and 
well understood - coagulation mechanisms are included, the process of grain growth is much too quick to be 
consistent with infrared observations of T Tauri disks. Small grains are removed so efficiently that, long before 
the disk reaches an age of 10^ years typical of T Tauri stars, the SED shows only very weak infrared excess. This 
is inconsistent with observed SEDs of most classical T Tauri stars. Small grains somehow need to be replenished, 
for instance by aggregate fragmentation through high-speed collisions. A very simplified calculation shows that 
when aggregate fragmentation is included, a quasi-stationary grain size distribution is obtained in which growth 
and fragmentation are in equilibrium. This quasi-stationary state may last 10*' years or even longer, dependent 
on the circumstances in the disk, and may bring the time scales into the right regime. If this is indeed the case, or 
if other processes are responsible for the replenishment of small grains, then the typical grain sizes inferred from 
infrared spectral features of T Tauri disks do not necessarily reflect the age of the system (small grains —> young, 
larger grains —> older), as is often proposed. Indeed, there is evidence reported in the literature that the typical 
inferred grain sizes do not correlate with the age of the star. Instead, it is more likely that the typical grain sizes 
found in T Tauri star (and Herbig Ae/Be star and Brown Dwarf) disks reflect the state of the disk in some more 
complicated way, e.g. the strength of the turbulence, the amount of dust mass transformed into planetesimals, the 
amount of gas lost via evaporation etc. A simple evolutionary scenario in which grains slowly grow from pristine 
O.l/im grains to larger grains over a period of a few Myr is most likely incorrect. 
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1. Introduction 

The coagulation of dust grains in protoplanetary disks is 
believed to be the first stage of planet formation. The dust 
grains, inherited from the interstellar medium, are initially 
particles of submicron size. But as a result of the high den- 
sities in the disk, and due to processes such as Brownian 
motion, settling and turbulence, the dust grains quickly 
get in contact with each other, and generally stick and 
form aggregates of ever increasing size. It is believed that 
the continuing growth of such grains/aggregates even- 
tually leads to the formation of planetesimals, which, 
through gravitational interaction, coalesce to form the 
rocky cores of planets. The subsequent accretion of gas 
onto these cores, if they are massive enough, then leads 
to the formation of gaseous giant planets. This is known 



as the core accretion model for giant planet formation 
fLissauer .199^) . 

A significant body of theoretical research has so far 
mostly focussed on our own solar system, trying to ex- 
plain the structure of the system as well as the rele- 
vant time scales known from geological and meteoritic re- 
search. Theoretical models are constrained for example 
by the fact that the core of Jupiter must have formed 
fast enough (within a few Myrs), to start accreting gas 
while this gas was still present in the disk, and to stop 
the Asteroid Belt from forming a planet. On the other 
hand, the spread in ages for some inclusions in mete- 
orites (CAIs and Chondrules) shows that that the pro- 
cess of putting t ogether planetesimals m ust have taken at 
least 4-6 Myrs ijBrearlev fc Jone "^. '1998*) . The pioneering 
work of Weidenschilhng l|l977n i980: 1984J considered the 
main processes contributing to particle growth, i.e. the 
collisons between grains caused by Brownian motion, ver- 
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tical settling, radial drift and turbulence. Many details 
in the process have since then been refined, such as the 
detailed physics o f the dust subla-yer for ming at the mid- 
plane of the disk toubr uUe et all llQQRj), the tr apping of 
particles in anti-cycloni c vortices (e.g. iKlahr fc H enning. 
Il997t ICuzzi et al.l[l993tl. the influence of aggregate shape 
on the timescales fWeidenschilling, '1997b*), as well as the- 
oretical (Dominik fc Tielens 1995, 1996, 199^ and exper- 
imental fe.g. lBhim & Wur'm[f200nl) studies of the strength 
and shape of the aggregates formed. 

One of the interesting results is a time-scale problem 
related to the radial drift of particles through the disk. 
Once a particle becomes of the size of cm or larger, it 
tends to decouple from the gas and move on Keplerian or- 
bits around the star. But the gas in the disk rotates with 
a slightly sub-Keplerian speed around the star, due to the 
small radial pressure support within the disk. The velocity 
difference between the gas and the particle causes friction 
and removes angular momentum from the dust particle. 
This results in a systematic drift inward toward the star. 
The time scale for this radial drift can be relatively short, 
on the order of 100 y ears for m-sizcd particles at 1 AU 
ifWeidenschillingj . Il977|) . Only once the particle has grown 
to kilometer size this radial drift stops because the friction 
has decreased sufficiently. The question is therefore: how 
can we grow grains from cm to km size within only 10^ 
orbits? Normal grain coagulation processes may or may 
not be fast enough, but particles of such large size have 
poor sticking properties. So the basic issue of grain coagu- 
lation is therefore how we can speed up coagulation in the 
models, and how can we enhance the sticking efficiency 
for larger particles. 

A completely different set of constraints on grain 
growth in protoplanetary disks can be derived from the ob- 
servations of circumstellar disks around young stars which 
are probably in the process of forming planets right now. 
The developments in infrared astronomy in recent years 
have opened up a new window to planet formation: the 
direct study of protoplanetary disks around T Tauri stars, 
Herbig Ae/Be stars and even around Brown Dwarfs. The 
ISO satellite, for instance, has provided infrared spec- 
tra of Herbig Ae/Be stars with unprecedente d acc uracy 
(e.g. Malfait et al. il999i: van den Ancker et al. 1200(f), giv- 
ing clues to the composition and size of dust particles in 
these disks (Bouwman et al. |2000) and to the geometry 
of these disks (Meeus et al. l2001il . The new Spitzer Space 
Telescope will do the same for T Tauri stars and Brown 
Dwarfs. With infrared interferometry (IOTA, PTI, VLTI, 
Keck etc) much has been, and will be learned about the 
dust miner alogy a nd structure of the disk at sub-AU scales 
(e.g. Eisner^S) 

and i n the ha bitable zones around these 
stars (van Boekel et al. l2004afl . All of these observations 
have drastically increased our knowledge of the physics 
of the birthplaces of planets. One of the main discover- 
ies from (sub-)mm observations (e.g. Testi et al. l2003|) 
an d from mid-infrared spectroscopy (e.g. van Boekel et 
is the ample evidence for grain growth in disks. 
However, the information about grain growth is encoded 



in a complicated way in the data we can obtain. The 
particles seen by different observations are different in 
size, and are located at different locations in the disk. 
Extracting the desired constraints on grain growth pro- 
cesses therefore requires self-consistent modeling of the 
growth processes with disk structure and radiative trans- 
fer. For the earliest phases (collapse and ear ly disk forma- 
tion) t his problem has been a ddressed by ISuttner et alJ 
(|l999l) : ISuttner fc Yorkj ((200lh . However, most observed 
star-|-disk systems have already long evolved beyond this 
very early phase. Also, the process of the formation of 
planets may span up to 10 Myrs or more, which is much 
longer than the duration of the early disk phases studied 
by Suttner et al. 

It is the purpose of this paper to confront self- 
consistent models of grain coagulation and settling in pro- 
toplanetary disks with the infrared observations of these 
objects. Interestingly, this has never been done before. In 
this paper we solve the coagulation equation (often re- 
ferred to as the Smoluchowski equation), coupled to the 
equation for grain settling and vertical turbulent mixing. 
The coagulation equation evolves the dust size distribu- 
tion in time, while the settling/mixing equation solves the 
vertical motion of the dust particles. At given time inter- 
vals we compute the spectral energy distribution (SED) 
and images of the disk using a 2-D axisymmetric contin- 
uum radiative transfer code. 

We will show in this paper that we are confronted with 
a new time scale problem, quite contrary to the one dis- 
cussed above: we find that the coagulation of small par- 
ticles is too fast to be consistent with infrared observa- 
tions. From infrared spectroscopy in the 10 micron regime 
it seems that grain growth from 0.1 /im to 2 /im hap- 
pens over a time scale of a few 10® years (van Boekel et 
al. ,2004b) . We will show with our models that it is very 
difficult to prevent the complete depletion of grains up 
to 100 /im within only lO'' years. We will suggest that 
aggregate fragmentation may provide a piece of the puz- 
zle. In doing so, we will necessarily repeat computations 
done before for the solar system, but we will do so with 
better resolution and with considering observational con- 
sequences. 

The structure of this paper is as follows: In Sect. El we 
describe the equations used to solve the problem of co- 
agulation in a protoplanetary disk. From there we build 
our model step by step, adding the different processes one 
at a time, in order to show the relative importance, and 
to demonstrate the solidness of the main conclusion. In 
Sect. 121 we first briefly discuss the fate of a single particle 
as it settles to the disk midplane and grows by sweeping up 
other grains. In Sect. 01 we solve the coagulation-settling 
equations for the entire ensemble of grains in a vertical 
slice in a disk, producing time-dependent grain size dis- 
tribution functions as a function of height above the mid- 
plane. In Sect. Owe use these local models to build a full 
disk model and to compute SEDs, with the conclusion that 
full coagulation depletes small grains much too rapidly in 
order to be consistent with infrared observations of disks. 
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In Sect. we introduce a possible remedy for this puzzle 
by showing that aggregate fragmentation might contin- 
uously replenish the small grains. Finally, in Sect. |7|we 
discuss our results. 

2. Equations 

The distribution of grains of a certain mass m, at a certain 
distance R from the star, a certain height Z above the mid- 
plane and at a time t is given by the distribution function 
f{R, Z, m, t). It is defined such that f{R, Z, to, t) m dm, is 
a dust mass density in g/cm^. The dust grains undergo 
several processes. First of all, they settle toward the mid- 
plane, and at the same time they get turbulently stirred 
back up again. In the absence of coagulation, an equilib- 
rium solution is eventually reached in which the upper 
layers of the disk are devoid of dust of a particular size, 
while below a certain Z the grains are fully mixed with 
the gas and hence have more or less constant abundance 
(e.g. DubruUe et al. 199fk Takeuchi & Lin 2002). This pro- 
cess of settling and stirring has pr ofound effects on the ap- 
pearan ce of protoplanetary disks l)Dullemond fc Dominikl . 
l2004bl) . 

As the settling and turbulent stirring process takes 
place, the grains also coagulate. The relative motions be- 
tween the dust grains, necessary for them to meet each 
other and form aggregates, are produced by various pro- 
cesses. In the models presented in this paper we include 
Brownian motion, differential settling and relative motions 
caused by turbulence. These are the most important pro- 
cesses for coagulation up to cm size particles, and we will 
describe in detail below how they are implemented in our 
model. 

The problem of dust settling and stirring is a problem 
in the coordinate Z while coagulation is a problem in the 
coordinate to. In principle they have to be solved simulta- 
neously. In practice, the settling and mixing are solved in 
the same subroutine, and the coagulation in another sub- 
routine. By using the technique of operator splitting one 
can switch between the two subroutines at each time step, 
thereby solving the simultaneous set of problems. 

2.1. Cross sections of particles 

An important quantity that will enter the equations is the 
cross section of a particle, for collisions with gas particles 
and with other particles. If we assume solid spheres, the 
collisional cross section for grain-gas collisions is simply 
given by CTg = 7ra^ where a is the radius of the particle. 
For collisions between two grains, the cross section will 
be cTc = 7r(ai -I- 02)'^. However, in reality, grains formed 
by aggregation are never compact spheres. The will have 
internal structure which depends on the formation mecha- 
nism. One extreme of the possibilities are Particle-Cluster 
Aggregates (PCA) which are formed when an aggregate 
grows by addition of small grains only. Such particles will 
tend to be spherical and porous, with a limiting porosity 
for large particles of about 90%. The other extreme are 



Cluster-Cluster- Aggregates (CCA) in which growth of an 
aggregate is dominated by accreting aggregates of its own 
size. Obviously, many other possibility exist, including si- 
multaneous mixture of PCA and CCA (growing by both 
small and large aggregates at the same time) and hier- 
archical growth by switching from one growth process to 
another, possibly several times. The purpose of this paper 
is not an in-depth study of these processes, but will be sub- 
ject of a future study. For the current study, we only use 
the three classical cases: compact, PCA and CCA parti- 
cles. For each of these growth classes, a unique relation be- 
tween mass and average collision cross sections can be de - 
rived. We here use the analytical fits bv lOssenkop3l)l993|) . 
These relations directly provide Og^irn) and <Tc(toi,TO2). 

2.2. Settling and vertical mixing 

The settling and mixing of grains is a 1-D time-dependent 
problem in z, which has to be solved for each R and to. We 
describe the settling and vertical mixing equa tions fo Uow- 
ing the discussion of Dullemond & Dominik ll2nn4bl) . We 
refer to that paper for details. The equilibrium settling 
velocity for particles smaller than about 1 cm (Epstein 
regime) is: 



Wsctt 
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where i^K = \/ GAI^/R^ is the Kepler rotational fre- 
quency, Cs is the isothermal sound speed. Ug is the col- 
lisional cross-section of the dust grain for collisions with 
gas or very small dust particles, i.e. the projected surface 
of the grain, averaged over all directions. 

The diffusion constant for turbulent mixing D is given 

by 



D 



(2) 



where aturb is the turbulence parameter. The symbol Sc 
is the Schmidt number defined as Sc = 1 -I- St with the 
Stokes number St given by 
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In this above equation g is a parameter characteriz- 
ing the turbulence. We take it q — 1/2 in this paper 
ISchraoler & Hcnning, 2004). 

With the vertical settling velocity and the turbulent 
mixing constant we can now define the settling/mixing 
equation as follows: 



dfim) 



dt 



dz 



d_ 
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In this paper we solve this equation for each radius R 
separately. We do not couple these radii, i.e. particles are 
not allowed to move from one vertical slice to another. 
In principle this should be included, since radial drift can 
be very important. But for simplicity, and since we are 
mainly interested in the growth from 0.1 /im grains to 1 
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cm grains, we ignore this process and treat every radius as 
a separate 1-D vertical settling/mixing problem, coupled 
to the coagulation equations described below. 

We solve the settling/mixing equation numerically 
using implicit differencing. This is necessary since the 
Courant-Friedrich-Lewy condition for the vertical mixing 
would require much too small time steps for an explicit so- 
lution, and would make the simulation prohibitively slow. 

2.3. Coagulation 

The coagulation of grains is a 1-D time-dependent problem 
in TO, which has to be solved at each R and Z. T he coagu- 
lation equation is (Schumann Eli^ TodesEiS Safronov 



df{m) 



dt 
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f{m')f{'m — m')fJc{m' , m — m!) 
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which is the continu ous fo rm of the Smoluchowski equa- 
tion (Smoluchowski Il916|) . The first term on the right 
hand side represents the gain of dust matter in the mass 
bin TO by coagulation of two grains of mass to' and to — to'. 
The second term represents the loss of dust matter in the 
mass bin to by coagulation of a particle of mass to with 
a particle of mass to'. Ai;(toi,TO2) denotes the average 
relative velocity between these two particles. This average 
velocity may consist of random motions but also of system- 
atic drift between particles of different mass. The combi- 
nation if(TOi,TO2) = (Tc(toi, TO2)Aw(toi, TO2) is Called the 
kernel of the coagulation equation. 

We include three processes leading to a Aw(toi,TO2): 
Brownian motion, differential settling and turbulence. For 
Brownian motion, the average relative velocity is given by: 



is very similar to the formation of rain drops in clouds in 
the earth atmosphere ('rain out'). The systematic relative 
velocity is: 



A'i;s(TOi, TO2) 
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This is clearly zero for particles with equal o^jm ratio. 
In the special case of a uniquely defined (Jg/m for any 
given TO, as we are using in this paper, this means that 
the relative velocity is zero for equal mass particles, since 
they both settle at equal speed. The differential settling 
as a source of collisions works best for particles of very 
different <Jg/m (and therefore mass). Therefore typically 
the largest mass aggregates will sweep up the smallest 
particles. This leads to relatively compact particles with 
porosities typically of the order of 90%. Such aggregates 
are often called particle-cluster aggregates (PCA). 

Turbulence-driven coagulation is a rather complex pro- 
cess. It requires detailed calculati ons of the statistics of 
motions of particles (Voelk et al. llSSOf) • Unfortunately, 
these calculations are too complex to directly build 
into a coagulat i on co de like the one we present here. 
IWeidenschillind l)l984|) has derived fitting formulae to the 
results of Voelk et al., which we implement with slight 
modifications here. We refer a discussion of these to ap- 
pendix^ It should be noted that these formulae introduce 
a considerable uncertainty into the present coagulation 
code, since they may depend on details of the turbulence 
which are not well known. For example, turbulence driven 
by th e magneto-rotational instability (Balbus & Hawlevl 
Il99l[) may have different properties than the turbulence 
assumed by Voelk et al.. The results we focus on in the 
current paper, which is the fast disappearance of small 
grains, is present even if we ignore turbulence as a source 
of relative velocities. Therefore, the uncertainties intro- 
duced by assuming a particular turbulence spectrum are 
acceptable. 

Solving Eq. ^ numerically is challenging. In appendix 
0we describe our numerical algorithm, and in appendix 
El we show the results of a comparison against a test case 
described in the literature. 



Aub(TOl,TO2) 



' 8fcT(TOi -|- TO2) 

7rTOiTO2 



(6) 



This represents an average of random velocities, which 
is highest when both particles have the smallest mass. 
Therefore, Brownian motion favors collisions in which at 
least one collision partner has low mass. Once these lowest 
mass particles get depleted, the next higher mass particles 
(which are aggregates of the smallest particles) will coag- 
ulate etc. Brownian motion will therefore lead to a nar- 
row size distribution which slowly moves to larger sizes. 
This hierarchical growth procedure leads to aggregates 
with a fracta l structure: so -called cluster-cluster aggre- 
gates (CCA) llMeakinLll987l) . 

Differential settling is the process by which large 
grains, which settle faster than small grains, sweep up the 
smaller grains on their way to the midplane. This process 



2.4. Grain fragmentation 

In most of this paper we shall ignore aggregate fragmenta- 
tion/destruction, since we are interested to see what hap- 
pens to the observables of the protoplanetary disk when 
coagulation is able to develop to its full extent. But in 
Section El we will present a single- vertical-slice calculation 
of a situation in which grain fragmentation is included. 
Including grain fragmentation in a proper way is challeng- 
ing. But since the purpose of that section is only to demon- 
strate the main effects of replenishment of the smaller 
grains, we are satisfied with an extremely simplified im- 
plementation of fragmentation. We assume two colliding 
aggregates to disintegrate entirely into monomers if their 
collision energy, when divided over the sum of the masses 
of the two particles, exceeds a certain value. As our frag- 
mentation energy we choose 1250 erg/g, which, for equal 
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mass particles, amounts to a collision velocity of 1 m/s. 
This represents rather loosely bound particles, which en- 
hances the effect of fragmentation. In reality somewhat 
higher fragmentation velocities may be required, and also 
other processes such as partial fragmentation (or crater- 
ing). But for the purpose of demonstration of principle 
this very simple recipe is justified. 

The recipe described above is implemented in the fol- 
lowing way. The mass bins below a < O.bfim are consid- 
ered to be 'monomer bins'. Any collision between parti- 
cles of size larger than 0.5/im with a collision energy (di- 
vided by mass) larger than the critical value will remove 
the mass from the two original bins and return this mass 
into the monomer bins. The distribution in which it is re- 
turned is fixed, and represents the starting distribution 
shape (MRN, see below). 

3. Prelude: the one-particle model 

Before we present our full-fledged coagulation-settling- 
mixing models we briefly revisit the sim ple on e-grain 
model discussed by Safronov in his book l|l969|) . since 
this model is very illustrative of the phenomena presented 
in the subsequent sections. The one-particle model fol- 
lows the growth and settling of a single dust particle in 
a disk, assuming that all other dust particles remain sus- 
pended in the disk and do not coagulate. As the particle 
settles, it sweeps up the small grains suspended in the 
disk. Therefore, the mass of the particle is an increasing 
function of time: ni{t) and the height of the particle above 
the midplane z(t) decreases. Here, and in the remainder of 
the paper, we take a very simple disk vertical structure in 
order not to complicate matters more than necessary. We 
assume that the disk is vertically isothermal and the gas 
density obeys hydrostatic equilibrium. We therefore have: 



■ exp 



( — 



where Hp is defined as 




(8) 



(9) 



with the mass of the central star, fi the mean molecu- 
lar weight, taken to be /i = 2.3 (for molecular gas) and R 
the distance from the star. The temperature T is assumed 
to be the midplane temperature of a disk irradiated un- 
der an angle of ai^r = 0.05 around a star of Af, = 0.5 Mq, 
T^, = 4000 K and i?, = 2.5 Rq. If we assume an isothermal 
disk structure without a warm surface layer, then the tem- 
perature of such a disk is T = a\l^^ rI / RT^ = 204 K. 

Now if we let the particle in question settle toward the 
midplane according to the settling velocity Eq. then 
the equation of motion for that particle becomes: 



dz 

—rr — ^sett 
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Fig. 1. The time-evolution of the height z, radius a and 
mass m of a dust grain in the simple one-particle model, 
for three different initial dust radii. The specific weight of 
the dust is 3.6 g/cm^. 



At the same time, the mass of the particle increases ac- 
cording to 



dm 
~dt 



p{z)\vsctt\uc{rn) , 



(11) 



(10) 



where the factor 0.01 is the dust-to-gas ratio. The above 
two equations form a coupled set of ordinary differential 
equations which can be integrated using standard integra- 
tion schemes. Since the evolving particle quickly is much 
larger than the particles suspended in the gas, we assume 
"■g = cTc - but we will see below and in section 01 that this 
is not fully correct. 

In Fig.nwe show the resulting z(f), m{t) and a{t) as- 
suming spherical compact silicate grains in a disk with 
S^gas = 10^ g/cvc? at i? = lAU, starting at a height 
z — A Hp (with Hp the pressure scale height of the disk 
given by Hp = ^ kT R? j \irapGM^ ) , and with grains of 
different initial radius (mass). As one can see from these 
figures, the grain grows exponentially as it sweeps up mat- 
ter during its decent. It reaches the midplane as a cm size 
pebble in a few hundred years, even though it would have 
taken a few million years to reach the midplane if the 
grain would not have grown to larger size on its way to 
the midplane. Interestingly the time it takes to reach the 
midplane is almost independent of the initial grain size. 
In fact, it is also almost independent on the compactness 
of the grains, as can be seen in Fig. |21 A porous grain has 
lower settling velocity but larger cross-section. It therefore 
sweeps up more matter, allowing it to regain the velocity 
it would have had when it would have been more compact 
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Fig. 2. As Fig. ^ but now for different porosities of the 
grains. 



l|WeidenschillineLll997b|) . Moreover, the porous grain ends 
up at the midplane with a larger mass than a compact 
grain. If one were to redo the experiment with fractal grain 
growth (resulting in cluster-cluster aggregates) the result- 
ing end size of the grain diverges completely, becoming 
(theoretically) larger than the entire disk but with virtu- 
ally infinite porosity. This is clearly unphysical. In reality, 
at some particle size a transition to more compact parti- 
cles with fractal dimension 3 will take place, avoiding this 
divergence. At what size and how this happens is still one 
of the main questions in the study of planet formation. 

As a preliminary conclusion from this simple exercise 
one can say that grain growth due to this 'raining effect' 
happens on a very short time scale, much shorter than the 
typical life time of protoplanetary disks. And it seems that 
porousness or fluffiness of the grain does not slow down 
this growth process. It merely increases the final grain 
mass. 

4. Local disk models with coagulation, settling 
and mixing 

The result of the previous section has shown that coagula- 
tion can proceed very quickly through the differential set- 
tling process. However, that calculation was based on the 
simple assumption that a single grain follows this mode 
of growth while all other grains remain suspended in the 
gas. In reality this is not the case. All grains evolve si- 
multaneously, and settle simultaneously. It is not clear 
at which height above the midplane this 'rain-out' pro- 
cess starts, and what the effect of the collective growth 
is. Moreover, as will be shown below, the initial mode of 
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Table 1. Table of parameters of single vertical slice 
models (the S-series). First column: model name, sec- 
ond column: Brownian motion, third column: differential 
settling ('rain effect'), fourth cohmm: turbulent mixing, 
fifth column: turbulence-driven coagulation and sixth col- 
umn: porosity (Compact, Particle-Cluster-Aggregate or 
Cluster-Cluster- Aggregate) . 



growth is Brownian motion, not differential settling, and 
only after a certain time the differential settling will take 
over. 

In this section we will show the results of the full- 
fledged coagulation-settling-mixing calculations following 
the equations outlined above. We do the calculation first 
for a single vertical slice at i? = lAU. As our disk vertical 
structure we adopt the same struct ure as in Sec t ion El 
Our initial distribution function is a lMathis et al.l I^1977^ 
(MRN) distribution from 0.1 /im to 0.5 /im. 

In order to show all the effects more clearly, we proceed 
in steps. First we assume no settling, nor mixing nor coag- 
ulation by differential settling nor coagulation by turbu- 
lence: we only include Brownian motion (model SI). Then 
we show a model in which we include the raining effect as 
well (the coagulation by differential settling), but still no 
vertical mixing (model S2). We then also include turbulent 
mixing (model S3). Finally we also include coagulation by 
turbulence (model S4) . These models will demonstrate the 
speed at which the coagulation takes place. 

The models SI, S2, S3 and S4 are for compact spher- 
ical silicate grains. In order to show what the effects of 
non-compactness would be, we also present models simi- 
lar to S4, but with particle-cluster aggregates (PCA) and 
cluster-cluster aggregates (CCA). These are the models 
S5 and S6 respectively. 

Tabic n provides an overview of the models of this sec- 
tion. The resulting midplane dust distributions at different 
times are are shown in Fig.OI In these plots, we show the 
actual mass distributions (i.e. firn) when plotted over 
logm or m • a • /(a) when plotted over log a, like vF^, for 
spectral energy distributions). On the x-axis we deliber- 
ately plot the radius of the grain a instead of the mass 
m, because the radius of grains is a more familiar quan- 
tity. For the PCA and CCA models (models S5 and S6) 
this radius is the equivalent radius as if the grain would 
have been compact (i.e. an equivalent radius a for the 
PCA/CCA models corresponds to the same particle mass 
as for the compact models). For the PCA/CCA models 
the real radius is much larger. 

For clarity, we plot the results in two different ways. 
Figure 01 shows the midplane size distributions at different 
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Fig. 3. The time-evolution of the distribution function for models SI- • • S6 (see tableQlfor the model definitions). 



times for the various models in a normal 2D plot. Figure01 
shows the same results in a 3D way, i.e. we provide a 
separate axis for time. This is useful in particular for the 
calculations in which the size distribution develops two 
different peaks - such development gets confusing in the 
2D plot. We plot m ■ a ■ f(a) as a function of logo so 
that surface under the distribution function corresponds 
to total density of dust at that position in the disk. 

As one can see from the upper left panels of Figs. 13 and 
01 the coagulation by Brownian motion tends to produce 
a peak distribution with a certain width. The peak moves 
toward larger sizes in a self-similar way. But the speed 
with which it moves toward larger sizes is proportional 
to apeak i^/^. This can be understood in very simple 
terms. If we assume for simplicity that we have to deal 
with a mono-disperse size distribution of compact parti- 



cles with size m(t) oc t^, it is clear that the number den- 
sity of particles is proportional to n{t) cx m 
The relative velocity is v{t) oc m cx For the 

compact particles in this calculation, the coUisional cross 
section is a{t) cx m{ty^^ cx t^f^/^. The change in mass of 
the particles is given by jit^'''^ = ^ cx m{t)n{t)a{t)v{t) . 
Comparing the powers to which t is raised in this equa- 
tion we easily find /i — 1 = /i — /x-|-2/i/3 — /i/2 with the 
solution /i — 6/5 implying m{t) cx t^^^ and a{t) cx t^/^. 
While for the very smallest grains this is an important 
growth mechanism, clearly this is not efficient to grow to 
very large sizes. 

If the differential settling is included (model S2) the 
initial stage of growth at the midplane is still Brownian 
motion. But at some point {t ~ 500 year) one can see a 
sudden 'rain shower' appearing around grain sizes of a ~ 1 
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mm. This differential settling has already started at higher 
elevations before that time, but has not yet reached the 
midplane until t ~ 500 year. The abrupt appearance of 
a ~ 1 mm size grains at t ~ 500 year means that the 
'rain drops' have finally reached the midplane and popu- 
late the midplane mass bins of size a ~ 1 mm, giving rise 
to the isolated peak distribution around this size seen in 
the upper-right panel of Fig. |2| During the 'rain shower' 
the height of the original peak of the distribution function 
(at smaller grain sizes) is strongly reduced: these small 
grains are swept up by the descending larger particles. 
This process is similar to what happens in the earth at- 
mosphere when aerosols (smog particles) are washed out 
of the sky by rain. 

The time scale for this 'rain shower' to reach the mid- 
plane is very similar to that predicted by the one-particle 
model. The particles are, however, about 3 times larger 
than in the one-particle model. This is because the col- 
lective raining allows grains to coagulate with particles of 
similar (though not equal) size, thereby increasing the col- 
lision cross-section of the grains. The maximum increase 
factor would be a factor of (Jc/cg — 4, if the grains were 
allowed to coagulate with equal size partners. But since 
the differential settling only works for particles of unequal 
size, this factor is reduced somewhat, i.e. becoming a fac- 
tor of 3 in our simulation. 

It is interesting to see that the peak value of to • a • / (a) 
of the rained-down grain population (after 500 years) is 
much higher than the peak value of the initial distribu- 
tion, even though the peak width is not much different 
from the width of the original distribution. This may ap- 
pear as a violation of mass conservation. The explanation 
is that due to the settling virtually all the dust grains 
larger than about 10 /im have settled in an extremely 
thin midplane layer with very high dust density. This is 
the "dust subdisk" in which, as is sometimes believed, 
gravitational instabilities may trigger planetcsimal forma- 
tion (Goldreich & Ward 1973; Youdin & Shu 2002). This 
thin midplane layer is allowed to form in model S2 be- 
cause we have explicitly switched off turbulence. In prac- 
tice there is doubt whether such a thin layer can exist, 
as even the slightest bit of tur bulence thickens the dust 
subdisk ijWeidenschillin 

After the 'rain-out' nothing much happens, since the 
'rain shower' has cleared out a large fraction of the small 
grains in the disk, and there is not enough of them left 
to continue the raining effect. Basically, the rain shower is 
over. The remaining population of small grains continues 
to grow via Brownian motion until it merges in the popu- 
lation of rained-down grains. It should be noted, however, 
that a minute amount of small grains from higher eleva- 
tions, which have survived the rain shower, still slowly but 
surely settles toward the midplane. These grains reach the 
midplane during the later stages of the simulations. The 
slightly larger grains arrive first, the smallest grains last. 
Once arrived at the midplane, the grains get incorporated 
into larger grains due to Brownian motion. The combined 
effect is a dip in the size distribution around about a size 



of 10/im. This dip is widening with time, producing a size 
distribution with two peaks. Due to the arrival of ever 
smaller grains, the small particle peak moves to the left, 
while the large particle peak slowly moves to the right 
(because of further Brownian motion grain growth) . If the 
small grains would not be incorporated into the big grains 
(for instance, if one would switch off Brownian motion co- 
agulation), the settled size distribution would freeze out, 
with the dip between the two peakes filled in. 

Around 1 Myr the secondary peak is seen at 1 /im. It 
appears as if this peak is only 10000 times smaller than 
the value of the initial distribution. While this is a strong 
depletion of small grains, it may not be enough to render 
the disk optically thin. However, here again this is some- 
what deceptive due to the thinness of the dust subdisk. 
For grains of this size the dust subdisk is 1% of the pres- 
sure scale height of the disk at 1 Myr. This brings the 
depletion factor to about 10^. 

If we include the vertical turbulent mixing with 
Q^turb = 0.01 (model S3), then a new interesting phe- 
nomenon occurs. The initial evolution remains the same as 
for S2 (Brownian motion followed by raining) , but the dif- 
ferential settling process does not stop anymore. Because 
the rained-down grains are now stirred up from the mid- 
plane again, they can rain down a second and third time, 
continuing to sweep up any grains they meet. This pro- 
longs the fast growth process considerably, and strongly 
depletes the small grains from the disk. This process is 
somewhat similar to the growth of giant hail stones in cu- 
mulonimbus clouds in the earth atmosphere: before such 
stones reach the earth surface, they often get transported 
back to high altitudes by the strong updraft within the 
cloud. In this way these hail stones get multiple chances to 
collect rain drops and smaller hail stones on their surface, 
allowing them to grow to sizes as large as decimeters. The 
main difference is that the air movement in cumulonimbus 
clouds constitutes a systematic flow, while in our models 
the turbulent motions are random. 

In model S4 we also include the coagulation by the tur- 
bulence itself. It is interesting to see that this process does 
not appear to make much difference in the earlier phases 
of the coagulation process. The differential setting effect 
(prolonged by the vertical stirring) is clearly the dominant 
process up to 10'* years. After that, the turbulence-driven 
coagulation takes over and manages to grow the grains 
even further. However, by this time we are well in the 
regime of boulders, in which many of the equations we 
use are no longer valid. Note a tiny wiggle in the distri- 
bution function at 10^ /xm: this is an effect that can be 
traced back to the discrete jump in the slope of the fitting 
formulae of appendix IXI 

From the above calculations it is clear that, if we al- 
low coagulation to work with perfect efficiency (i.e. maxi- 
mum sticking, no destruction) coagulation happens on an 
extremely short time scale, even if we only include coag- 
ulation by Brownian motion and differential settling into 
our calculations. If we would include other processes such 
as coagulation by radial drift, while still assuming perfect 
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Fig. 5. The vertical optical depth in the UV at 1 AU (i.e. 
the integrated projected surface of the dust particles) as 
a function of time for the vertical slice models SI. . . S6. 



efficiency, then this is only aggravated instead of allevi- 
ated. It seems therefore that this result is rather robust. 
We shall see in section El that destructive collisions may 
indeed play an important role. 



4.1. The optical depth of the models 

A first hint on the effects of the coagulation on the ap- 
pearance of a disk can be taken from the optical depth at 
short wavelength, i.e. at wavelength corresponding to the 
stellar photons which are heating the disk surface and are 
responsible for the flaring in the disk. 

In Fig. 13 we show the vertical optical depth at UV 
wavelengths through the disk. For the cross section we 
have taken the projected geometrical cross section, i.e. 
(Tg. We can see that for all models, the optical depth de- 
creases significantly due to the coagulation. While in the 
pure Brownian Motion case SI, the optical depth decreases 
by a factor of 100 in lO^yrs, including settling already 
decreases the optical depth by more than four orders of 
magnitude in the same time, pushing the optical depth be- 
low 1. With turbulent mixing and turbulent coagulation, 
the effect becomes enormous, and the disk becomes fully 
transparent in only lOOOyrs. Using PCA and CCA grains 
does slow down the process somewhat. In particular, in 
the CCA calculation the optical depth stays constant for 
about 1000 years, but then also starts to decrease quickly. 
After 10^ years, also in this case the optical depth is be- 
low one. The possibility of change i n optical depth o n short 
time s cales was also indicated bv IWeidenschillinS l)l98Ct 
Il997al) who found that the Rosseland optical depth may 
decrease by an order of magnitude in about 1000 orbital 
periods. 



Table 2. Same as tabled but now for the full disk models 
Fl and F2. 



5. Global disk models with coagulation, settling 
and mixing 

The short coagulation time scales will have consequences 
for the infrared and optical appearance of protoplanetary 
disks. In this section we glue a series of vertical slices (an- 
nuli) together to form a full disk model. We perform the 
coagulation-settling-mixing calculations in each of these 
slices and make snapshots of the distribution function 
at given times. From this series of slices we construct 
a 3-D axisymmetric (i.e. effectively 2-D) disk model at 
each of the snapshot times, and feed these disk mod- 
els into a multi-dimensional continuum radiativ e trans - 
fer code called RADMC (see Dullemond & Dominik l2004a|) . 
With this code we can produce synthetic spectral energy 
distributions (SEDs) for the disk at the given times, as 
well as images at various wavelengths and inclinations. 

Our model has an inner radius Rm = 0.7 AU, outer 
radius i?out — 100 AU and has a gas surface density 
profile i;gas(i?) = So(i?/AU)'' with p = -1.5 and 
So = 400g/cTO^, which amounts to a disk with Af^isk — 
0.005 M0. In total we glue 40 vertical slices together to 
form the full disk model. The radii of these slices are log- 
spaced between Rm and i?out, meaning that the width of 
each slice is about 14% of its radius. We take the same 
stellar parameters as in previous sections, and similarly 
we assume compact silicate grains. The temperature of 
the disk is a function of radius R only, i.e. it is constant in 
vertical direction. As in the previous sections, the temper- 
ature is assumed to follow from thermal balance assuming 
that the disk is passive and irradiated by the star under an 
angle = 0.05, leading to T{R) = a\l^^ ^ R^/ RT,.. We 
present two models. Model Fl has no turbulent stirring 
(we set aturb = 0) while model F2 has turbulent stirring 
with aturb = 0.01 (see table EJ. 

The left column of Fig. El shows the evolution of the 
SED of the full disk model without turbulence (Fl). One 
sees that very early the mid-IR drops, while the far-IR re- 
mains. This is because the coagulation processes are faster 
at small radii, and therefore the depletion of small grains 
happens faster at small radii than in the outer regions of 
the disk. It is also clear that the 10 /xm feature quickly 
loses strength but does not disappear before most of the 
mid-IR continuum flux disappears as well. Its shape also 
does not change appreciably. This is an interesting phe- 
nomenon, since observations of the 10 /zm feature of T 
Tauri stars and Herbig Ae/Be stars often show flatten- 
ing and weakening of this feature which is generally in- 
terpreted as a signature of gra in growth (v an Bo ekel et 
al. 2003; Przygodda et al. l2003t Meeus et al. l2003t Honda 
et al. 1200311 . According to the present calculations such 
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Fig. 6. The SEDs (left column and surface heights (right column) of the full disk models. The two top panels show 
the result for model Fl, the calculation without turbulence. The bottom panel results are for model F2, including 
turbulence. Different lines indicate the results after different times, as shown in the legends. 



flattened features are not predicted. Moreover, the entire 
mid-IR flux vanishes much too quickly (well within 10^ 
years). In fact, by 10® years most of the IR excess (also 
far-IR) has vanished, which is clearly inconsistent with ob- 
servations of T Tauri stars and Herbig Ae/Be stars. The 
same behavior is reflected by the surface height in the disk, 
plotted in the upper right panel of figure El The surface 
height is calculated by following starlight radially away 
from the stellar surface and determining at what location 
an optical depth of unity for a wavelength of 0.55 fim is 
reached. In contrast to the vertical optical depth shown 
for the slab calculations (see figure ]^ , this plot is also 
sensitive to the vertical distribution of the material. It 
measures both loss in total optical depth, and settling of 
the opacity carriers toward the midplane and therefore is 
the most important indicator for understanding the SED. 
The surface height initially begins to decrease globally, 
which is mainly due to the settling motion of particles in 
the uppermost layers jlSullemond fc Do minik. 2004b) . In 
the inner disk regions, the effects of coagulation quickly 
lead to a transparent disk, with the disk "surface" dis- 
appearing. The boundary where photons are intercepted 



moves outward and reaches 30 AU in 10® years, while in 
the outer disk, the surface height continuously decreases. 
If one includes turbulence in the model, then the SEDs 
become as shown in the two bottom panels of Fig|3 In 
this case the problem is even more acute. 

It is interesting to see that as the mid-IR flux van- 
ishes, the far IR flux temporarily increases. This is be- 
cause of energy conservation. The dust coagulation hap- 
pens much faster at small radii than at large radii. As the 
inner regions of the disk become optically thin, the outer 
regions are still optically thick and reprocess all the radi- 
ation that was previously reprocessed by both the inner 
and the outer regions of the disk. Once also the outer re- 
gions get deprived of their small grains, they too become 
optically thin and the entire IR excess drops down. Note 
that the far-IR excess for model F2 (with turbulence) is 
usually larger than that of model Fl. This is because even 
though the grains grow faster in model F2, they also get 
swept up higher above the midplane so that the disk can 
reprocess a larger fraction of the stellar radiation. 

The results of this section clearly show that the quick 
disappearance of the small grains due to the efficient inter- 
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where it immediately starts to coagulate again. For grain 
sizes larger than 1 cm the distribution function continues 
to evolve, forming a powerlaw 'tail' with a growing peak 
value. This is a tail distribution of 'lucky' grains which 
managed to avoid an encounter with an equal size colli- 
sion partner and therefore avoided fragmentation. These 
grains sweep up some of the small grains from the semi- 
stationary distribution below 1 cm, and therefore manage 
to grow toward larger sizes. By this time other important 
effects, which we haven't taken into account yet, will be- 
come dominant, such as radial drift and the consequent 
run-away growth. Also some of our equations start to be- 
come invalid for such large particles, in particular our for- 
mula for the gas drag in the Epstein regime. 



Fig. 7. The time-evolution of the distribution function for 
model SDl (including all growth processes and aggregate 
fragmentation). This sUce has only been evolved up to 
8 X 10^ years (darkest curve). 

action between the differential settling and the turbulent 
mixing is clearly in contradiction with observations. The 
absence of turbulence may leave a minute population of 
small grains, but is not very efficient in solving the prob- 
lem fully. In the next section we discuss what we believe 
is the most likely solution to this apparent contradiction 
with observations. In Section [7| we discuss various other 
possible solutions. 

6. A simple local model including aggregate 
fragmentation 

So far we have taken a perfect sticking condition: if two 
grains/aggregates collide, they stick and form a larger ag- 
gregate. However, for large collisions velocities this is no 
longer true. In this section we show what happens when 
aggregates are allowed to disintegrate upon collision. This 
replenishes the small grain size bins, and may go some way 
toward solving the time-scale problem. Unfortunately the 
computational demand for such a calculation is orders of 
magnitude higher than for a pure coagulation calculation. 
Therefore we only show the results for a single vertical 
slice (model SDl). 

We include aggregate fragmentation in the very sim- 
plified way described in Section if the collision energy 
exceeds a certain limit, both grains are destroyed upon 
impact. We take as our model parameters again the same 
parameters as in model S4 of Section^ 

From Fig. |7| we see that, while the initial stages of 
growth are identical to those of model S4, very quickly the 
replenishment of small grains due to fragmentation starts 
to take place. After about 10'' years a semi-stationary state 
is reached for sizes below a < 1cm. This is an equilibrium 
between grain growth and grain fragmentation. As grain 
aggregates reach sizes of about 1 cm, most of the aggre- 
gates are destroyed and their mass returned to monomers. 



7. Discussion 

7.1. Other mechanisms to keep up the small particle 
population 

The basic fact that coagulation will lead t o a reduced opti- 
cal depth in disks has been noted earlier. IWeidenschillind 
(^984) noted that the reduction in opacity may lead to 
a termination of turbulent gas motions i f thos e motions 
are driven by convection. iMiz^mo et al.l ((l98? ^ forced a 
steady-state solution for their coagulation calculations by 
assuming that the disk is constantly supplied with new 
gas containing new small particles. They also discuss that 
the addition of new particles may lead to alternating 
phases of high and low optical depth, and consequently of 
convection driven turbulence. However, while turbulence 
driven by convection may be present in disks, other drivers 
for turbulence may b e important as wel l . The magneto- 
rotational instability ijBalbus fc Hawlevl Il99l|) seems to 
be an important candidate. For the presence of this in- 
stability, low gas column densities are required in order 
to allow cosmic rays and X rays to penetrate and ionize 
t he gas. High ga s columns can create so-called dead-zones 
Gammielll996|) in the disk mid-plane where ionizing radi- 
ation (such as X-rays and cosmic rays) cannot penetrate. 
Furthermore, the idea of a constant inflow of small dust 
grains onto the disk can only be valid for the early evolu- 
tionary phases of a disk. Observations show that circum- 
stellar disks at a n age of several milli on years still can be 
strongly flaring l)Leinert et all l2004|) . By this time, the 
parental cloud has been largely cleared away and infall of 
new material as well as accretion onto the star has vir- 
tually ceased. The exact limit on the amount of infalling 
material which can be relevant for the disk opacity is not 
easily estimated. It will depend on the speed at which ver- 
tical mixing and turbulence can remove these grains from 
the disk. Clearly, for newly added material, the removal 
time scale will be smaller since the low dust density slows 
down coagulation. Pending more detailed calculations of 
these effects, it seems to us that the infall of fresh gas and 
dust will not solve the discrepancy between the speed and 
effectiveness of coagulation on the one hand, and infrared 
observations of disks on the other hand. 
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7.2. The effect of the large grain population on the 
small grains 

The choice of upper boundary to the grain size (omax) in 
the above presented simulations has a strong influence on 
the results for the case in which both differential settling 
and turbulent mixing are included into the calculation. 
If flmax is taken too small (for instance 1 mm only), the 
grain growth is artificially stopped at that size. This means 
that a comparatively large population of grains remains at 
a size around Omax (while it should have grown further), 
providing a comparatively large cross-section for the de- 
pletion of the smaller (1 fj,m) grains. Therefore, choosing 
flmax to be 1 mm instead of 10 m or more, strongly en- 
hances the depletion of the small grains. In our simula- 
tions we therefore had to choose amax to be very large 
(we took it 100 m). But this introduces yet another prob- 
lem: our equations are no longer accurate at those large 
sizes. Moreover, large bodies tend to drift inward toward 
the star at a very high pace. This is not included in the 
present model, because each vertical slice is assumed to be 
independent of the other slices. Future investigations will 
have to deal with this problem in a proper way. However, 
it is not likely that that would solve the problem of the 
quick depletion of small grains, since even for the non- 
turbulent case (which represents the slowest growth) the 
depletion is quite strong. 

7.3. The inner regions of the disk 

Close to the central star (the inner regions of the disk) 
the depletion of small grains goes the quickest, since the 
Kepler time scale is the shortest there. The problem of 
the quick depletion of small grains seems therefore to be 
most acute in those regions. But can we be sure that small 
grains indeed exist so close to the star? The answer seems 
to come from very recent observations with the mid in- 
frared interferometer MIDI on the Very Large Telescope. 
Using this instrument, van Boekel et al. (Nature, submit- 
ted) separated the correlated 8-13 /zm flux from the total 
8-13 iiui spectrum for 3 Herbig Ae/Be stars. In this way 
they were able to single out the spectrum from the in- 
ner 2 AU region of the disk. Although the 10 fim silicate 
feature in this correlated flux spectrum clearly showed ev- 
idence for grain growth up to 2 /j,m, it also clearly showed 
that grains of approximately 2 /im were still present and 
their abundance is strong enough to produce a clearly 
discernible 10 /xm feature. This is in clear contradiction 
with the pure-coagulation models presented here, and re- 
inforce our conclusion that aggregate fragmentation (or 
some other resupply of small grains) should play a major 
role in disks. 

7.4. Interpretation of the results 

The models of grain growth presented in this paper show 
that coagulation happens on a time scale that is two to 
three orders of magnitude too short to be consistent with 



observations of T Tauri star disks or the disks around 
Herbig Ae/Be stars. Turbulent mixing combined with co- 
agulation through differential settling is highly efficient 
at removing the small grains from the disk at all heights 
above the midplane. One may argue that some disks may 
have regions of zero turbule nce (f or instance the 'dead 
zone' introduced by Gammie '1996') , which may decrease 
the efficiency of the grain growth. But even if we only have 
differential settling, but no vertical mixing (as in model 
S2) then only a tiny population of small grains remains 
in the disk, containing less than 10^^ of the original pop- 
ulation of small grains. The full disk model Fl, which is 
without turbulence, clearly shows that at typical ages of 
T Tauri stars (around 1 Myr), almost no IR emission is 
left, and for the model F2 (with turbulence) this is even 
more dramatic. If anything, the SED looks like that of a 
debris disk instead of a T Tauri star disk. But even if we 
compare the SED at 10"* years to observed SEDs, we find 
that the dip in the SED at near- to mid-IR wavelength is 
rather untypical for most T Tauri stars. 

There are, however, a few examples of objects which 
show conspicuous near/mid-IR dips. A strong near/mid- 
IR dip in the SED of TW Hydra was interpreted by Calvet 
et al. l)2002|) as a signature of a planet which has cleared 
out the inner disk. The synthetic SED of model F2 at 10** 
years, however, shows a similar dip, and therefore dust co- 
agulation could be an alternative explanation. HD100546 
is a Herbig star which also features a conspicuously weak 
near-IR excess, which was a ttributed to a huge gap in the 
disk (Bouwman et al. l2003l) . Here again, dust coagulation 
may be an alternative explanation. 

From the results presented in this paper it seems un- 
avoidable that some form of replenishment of small grains 
is needed to make the model calculations comply with the 
observations. The only other possibility is that the stick- 
ing probability is enormously reduced by some process. 
Since we are not aware of a process capable of reducing the 
sticking probability by such a dramatic amount, we believe 
that replenishment is the only solution. Replenishment by 
destructive collisions seems to be the most natural way 
to prevent the small grains from disappearing entirely. In 
this paper we demonstrated that this could work if we as- 
sume very low binding energies of the grains. The process 
of cratering (a small particle impacting on a bigger one 
and creating a certain amount of impact debris) may be a 
better way to produce small grains, but we defer a more 
detailed implementation of aggregate fragmentation to a 
future paper. 

8. Conclusions 

In this paper we have modeled dust coagulation in proto- 
planetary disks and computed the SEDs of such disks. We 
qualitatively compare these SEDs to what is typically ob- 
served from T Tauri star and Herbig Ae/Be star disks, and 
we conclude that if coagulation is allowed to proceed un- 
hindered (i.e. without fragmentation of aggregates), then 
the small grains are depleted on a time scale that is three 
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orders of magnitude too short to be consistent with these 
observations. We have included three coagulation mech- 
anisms in this model (Brownian motion, differential set- 
tling and turbulence). The inclusion of only the first two, 
well understood, processes already shows that the strong 
and rapid depletion of small grains is unavoidable unless 
small grains are somehow replenished. The inclusion of 
a little bit of turbulent mixing will only aggravate mat- 
ters. It is very difficult to slow this process down, even 
by grain charging or other mechanisms. Either the grain 
sticking efficiency is many orders of magnitude less than 
currently assumed, or the small grain population must 
bo replenished in some way. We suggest that aggregate 
fragmentation could be such a mechanism. We present a 
highly simplified model for this and show that a semi- 
stationary equilibrium sets in in which coagulation and 
fragmentation are balanced for an extended amount of 
time. Whether this is the solution to the paradox remains 
unclear and requires much more detailed simulations. 
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Appendix A: Turbulence-driven coagulation 

Turbulent motions in the gas can cause collisions between 
dust particles. The basic reason for this is that particles 
with different ig/m ratio couple to eddies of different size 
and therefore acquire random velocities. To calculate the 
motions caused by turbulence, it is necessary to integrate 
over the contributions of all eddies from the largest scales 
down to the smallest scales which are set by the condition 
that the turbulent Renolds number Re equals 1. 

Calculating the gas-dust interaction is involved 
and has been covered by previous authors (e.g. 
IWeidenschiiih^ Il977t ICuzzi et all Il993() . Here we only 
describe the basic recipe implemented in our code. 

The main particle property that enters the calculation 
is the stopping time of of a particle which is given by 



3ml 
4 CTg pCs 



(A.l) 



which physically is the time in which a particle reacts to 
changes in the motion of the surrounding gas. Specifically, 
for a given turbulent eddy, this time scale indicates if the 
particle will follow the eddy motion or not. 

In order to derive the average relative velocities re- 
sulting from this mechanism, the full structure and spec- 
trum of the turbulence must be known. Turbulence is often 
characteriz ed by the velocity and i: i me sca le of the largest 
eddies (see IPuUemond fc Dominikl . l2004bl) 
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where g is a turbulence paramet er be tween and 1. 
Following Schraepler and Henning l|2004fl we take it to be 
q = 1/2. The energy then cascades down to small sizes, 
until the flow becomes laminar at a turbulent Reynolds 
number Re= 1 in a gas with viscosity 
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(1980*) and iMizuno et aP l|l988|) . IWeidenschillh^ I^1984^ 

has fitted the numerical results with a simple analyti- 
cal formula, whi c h has also been used for example by 
ISuttner k Yorkd l)200l|) . We wiU also adopt it, but with 
two modifications: 



1, 



Because of the lower cutoff to the eddy size spectrum, 
no random velocities relative to the gas are introduced 
by turbulence for particles with a stopping time below 
the turnover ti me of the smallest eddy , t^^^j. This has 
been noted bv IWeidenschilfind I^1984^ and discussed 
in detail bv iMizunoTr^iLl l)l988|) . We use the limiting 
case provided by the latter authors. 
The analytical fit produces an overshoot in the 



limit t2 



If U <C to, the analytical fit 



given by IWeidenschillind l)l984|) and also used by 
ISuttner fc Yorkd l)200l|rieads to v°^. = 3z ;L^ while 
the numerical results by IVoelk et al.l l) l980^ only ex- 



ceed 



^cdd 



by a few percent. In order to avoid effects 



cause by this incorrect high collision speeds, we there- 
fore limit the relative velocity caused by turbulence to 



odd- 



With both modifications, the recipe for the turbulent 
collision velocities becomes 
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Appendix B: Numerically solving the coagulation 
equation 

The numerical solution of Eq. |SJ| is a subtle matter. 
Consider a discretized grain mass grid rrii with i g [l^N]. 
The first integral on the right hand side can be con- 
verted into a sum over m' = (for fixed m = rrii) with 
k = \ . . .1 where / is the highest index for which mi < rm. 
Unfortunately the value of rui — ruk lies generally not ex- 
actly at a discrete mass point and therefore the value of 
/(m — m') = f{mi — nik) must be obtained by interpo- 
lation. Since / can vary extremely strongly, the interpo- 
lation is best done in log(/) instead of in /. Numerical 
practice has shown that two-point interpolation makes the 
algorithm more stable than four-point interpolation. The 
last point in the integral (i.e. the numerical sum) is lo- 
cated at m/2, and in that case ni' = m — m', i.e. both m' 
and m — m' are located in between two mass points and 
similarly an interpolation needs to be done. The integral 
in the second term on the right hand side of Eq. ^ is 
easier, since no interpolation is needed. 



The random motions of the particles must then be cal- 
culated by integrating over the contributions of the differ- 
ent eddy scales. For a Kolmogorov spectru m of the tur- 
bulence, this has been done numerically by IVoelk et al.l 



B.l. Renormalization 

The right-hand-side of Eq. ^ consists of a gain and a loss 
term. The gain term describes how much matter enters a 
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certain mass bin through coagulation of smaller particles, 
while the loss term describes how much matter leaves the 
mass bin through coagulation of this matter with particles 
of any other size. Typically these two terms are very large 
numbers which cancel each other almost entirely, except 
for a tiny amount. It is this tiny amount that is the cru- 
cial source term for the coagulation equation. This near 
cancellation happens then when the gain of matter in a 
bin is dominated by coagulation between large and small 
particles. The reason for this near cancellation is best de- 
scribed with an example. If a rock of 1 kg that hits a 
dust particle of 1 micron size, formally the rock increases 
in mass (albeit by an extremely small amount). The rock 
is therefore removed from its mass bin and put into the 
mass bin a few picogram toward larger mass. Since the 1 
kg rock may collide with trillions of 1 micron size parti- 
cles, the gain and loss terms in the kg mass bin are huge, 
but virtually identical. The minuscule difference between 
these two terms determines the eventual growth from 1 kg 
to 2 kg after the rock collects 1 kg worth of micron size 
particles. 

Numerically this poses a significant challenge. If one 
simply computes the gain and loss terms, the near cancel- 
lation goes astray once the cancellation happens beyond 
the 14th digit. In effect the near cancellation turns into a 
perfect cancellation, which is incorrect. To solve the prob- 
lem the integrals of the gain and loss terms have to be cal- 
culated simultaneously. The integrands of both integrals 
are calculated at the same time, and then subtracted and 
collected into a single integral. At each m' it is checked 
if the two terms tend to produce a near cancellation. If 
so, their difference is recomputed using a renormalization 
technique: 



/(m')/(rn — m')(Jc{m' , m — m')Av(m' , m — m') — 
f{m')f{m)ac{'m'^ r7i)Ati(m', m)dm' 
~ — f{m!){m, — to') 

(B.l) 

d(/(TO")ac(TO', m")Az;(m', m")) 



dm" 



which uses rHopital's rule. This renormalized version of 
the integrand is valid only when a near-cancellation takes 
place, and remains numerically well-determined even for 
extreme cases of near-cancellation. 



B.2. Mass conservation 

By performing the integration in the above described way, 
total dust mass is not necessarily perfectly conserved. 
Small errors can increase or decrease the total mass, and 
since these errors are generally systematic, the risk is high 
that the dust mass will unlimitedly grow or diminish as 
the simulation progresses. In our code this is avoided by 
making small corrections. Denote S{m) as the right-hand- 



side of Eq. ©. Define now the following two functions: 

, X _ / where S{m) > 
O-^TOj - \s{m) where S'(m) < ^ ^ 

The new right-hand-side S'ncw("^) is now: 



S'now(TO) = S-{m,) + xS+{m) 
where x is defined as: 

/ TO S- (to) dm 



X = 



J TO 5*4. (to) dm 
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This is generally a tiny correction, and it guarantees mass 
conservation and thereby prevents unphysical build-up or 
loss of matter. 



B.3. Time step determination for the coagulation 

Since the evolution of the grain size distribution proceeds 
each time step by adding a source term with the results 
of the integrals, the magnitude of this time step is set by: 



At = ^ min I ^ 



\\S{m)\J 

where < ^ < 1 is an accuracy parameter which we usu- 
ally set to 0.3 in our simulations. Typically this time step is 
shortest at the smallest to, so that the evolution at small to 
limits the time step of the entire simulation. This could in 
principle be prohibitively short compared to the total time 
we wish to evolve our simulation (10'' years). Fortunately, 
as the aggregates grow, the smallest mass bins are quickly 
depleted and can be set to zero once the value of /(to) 
drops below some floor value. The time step then does 
not need to be limited anymore by these small mass bins. 
In practice this means that as the time progresses, the 
time step becomes larger and 10^ years can be reached 
without many problems. 

Sometimes, the simulation may get temporarily stuck 
at relatively small time steps, despite the above men- 
tioned time stepping method. This is caused by the op- 
erator splitting between the settling/mixing and the co- 
agulation. While the coagulation equation may attempt 
to empty a certain mass bin, the settling and mixing may 
fill it up again. This typically happens at the midplane, 
where settling and vertical mixing may replenish the mass 
in a certain mass bin by transporting it down to the mid- 
plane from higher altitudes. Since the settling and mixing 
are numerically simulated in an implicit way, allowing the 
time step to exceed the Courant condition for for settling 
and mixing by many orders of magnitude, this can cause 
the instant refilling of the mass bin after the coagulation 
equation has tried to deplete it. In practice, however, the 
simulation never gets entirely stuck, and at some point 
manages again to increase the time step to large enough 
values to allow the simulation to end in about 20 minutes 
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per vertical slice on a Pentium 4 processor. Therefore it 
is manageable to simulate an entire disk consisting of 20 
vertical slices in about 8 hours CPU time. 

However, once the aggregate fragmentation is in- 
cluded, the small grains remain present, and the time step 
remains limited by the smallest mass bin. The above time 
stepping method then does not work anymore, and the 
simulation may take large amounts of CPU time even for 
a single vertical slice. A fully implicit treatment of the 
coagulation/fragmentation may then be necessary to pre- 
vent excessive computational costs. 

Appendix C: Test cases 

It is not straightforward to test the numerical algorithm 
described here with the physics we include into the coag- 
ulation equation, since n o anal ytic solutions exist for the 
kernel we use. Lai et al. have presented analytical 

solutions to the coagulation equation with a Brownian mo- 
tion kernel, but these solutions are only valid in limiting 
cases and they contain undetermined constants. 

For simplified kernels, various analytic solutions to 
the coagulation equation exist. The two most well-known 
cases are the case of K{mi, 7712) — {mi + 7712) A (Safronov 
Il969t henceforth test 1) and the case of K{mi,m2) — A 
(Smoluchowski 1916; henceforth test 2). In both cases A 
is an arbitrary constant which we take to be A = 1 . Bot h 
tests were extensively discussed by Ohtsuki et al. l|l990|) . 
They find that in particular for test 1 numerical algo- 
rithms generally deviate significantly from the analytic 
solution if the coordinate spacing in mass is too coarse 
{rrii+i/mi > \/2). This seems to be a particularly tough 
test case. In contrast they find that test 2 is much less 
challenging for numerical algorithms. 

We have done both tests with our algorithm, follow- 
ing Ohtsuki's test procedure and using the analytic solu- 
tions presented in that paper. We find similar results as 
they find. For test 1 we get qualitative agreement, but 
the precise location of the peak of the distribution seems 
to depend on mass grid resolution and time step size. In 
Fig. ED the results are shown at a single given time, for 
four simulations: for 100 and 1000 grid points and for the 
largest possible time step (Aimax = 0.3max(/(m)/5'(?7i)), 
where the factor 0.3 is our "safety parameter") and for 0.1 
time that value. The modest but non-negligible deviations 
found in test 1 are slightly troubling. But since Ohtsuki 
et al. experience similar problems and the kernel of test 
1 is very unlike the kernel we use in our models, we are 
reasonably satisfied with the results of these tests. 

For test 2 we get excellent agreement for moderate 
grid resolution (100 grid points) and for the largest possi- 
ble time step size described above. The results are shown 
in Fig. lC.2l Initially the difference between our model and 
the analytic solution is relatively large, but this is be- 
cause our coagulation equation is formulated in a contin- 
uous form while the Smoluchowski solution is based on 
the discrete form of the equation. We therefore start with 
a different initial condition as in the analytic solution of 




Fig. C.l. A snapshot of the distribution function of test 1 
at dimensionless time r = 9, for various values of the grid 
resolution and time step. The symbol At is the time step 
{Atniax = max(/(77i)/5'(m))) and N is the number of grid 
points in m {N = 100 corresponds to mi^i/rrii — \/2, 
i.e. the highest resolution used by Ohtsuki; N = 1000 
corresponds to mi+i/rrii — 1.03517). The solid line is the 
analytic solution of Safronov. 
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Test 2 (Smoluchowski Test) 
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Fig. C.2. The time evolution of the distribution function 
in test 2 for the model (solid line) compared to the analytic 
solution of Smoluchowski (dotted line) . The snapshots are 
at dimensionless time r = 1, 10, 10^ 10^, 10", 10^ from left 
to right. Due to different initial conditions the agreement 
is not good at r = 1 but becomes better as time progresses 
because the initial conditions are 'forgotten' by the sys- 
tem. The tiny wiggle in the solid curve around to = 4 is 
a remainder of the initial conditions that apparently the 
algorithm does not 'forget'. 



Smoluchowski and logically we get initially somewhat dif- 
ferent results. But as time progresses, the initial conditions 
are 'forgotten', and the model result and the analytic so- 
lution start to agree better and better, ending with almost 
perfect agreement in our last time step. 
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These positive test results for simple kernels give hope 
that the full code, with dust settling, mixing and coagu- 
lation through Brownian motion and differential settling, 
also works well. Moreover, the stability and reliability of 
the code was confirmed by reproducing the results pre- 
sented in this paper with different grid resolutions (both 
in z and in m), and by using different time steps. But of 
course the only way to be sure that the code is working 
properly is a comparison of our code with various inde- 
pendent codes written by independent authors. Such a 
comparison has not been carried out. There is, however, 
an interesting way to test the code independently: by com- 
paring the average size of grains rained out onto the mid- 
plane (in the absence of turbulence) to the result of the 
one-particle model. From the one-particle model (Section 
it follows that the size of the grain, once it reaches the 
midplane, is a = 0.01 S/Sp^ for compact grains. It should 
be kept in mind that this result was obtained by assuming 
that only the test particle rains down, while the other par- 
ticles remain suspended in the disk. We can simulate this 
effect in the full code by taking two populations of grains: 
a fixed population (equal to the initial grain population) 
and an evolving population. The evolving population only 
grows by colliding with the fixed population, while the 
fixed population is not changing during the simulation. 
This does not conserve mass, but it does simulate the con- 
ditions similar to the one-particle model. We carried out 
this test for Egas = 100 g/cm^ at i? = lAU with pa = 3.6, 
and we find that the model produces a sharp peak dis- 
tribution near a = 0.35mm, which is virtually identical 
to the analytic result (a = 0.3472mm). That it is slightly 
larger is due the fact that the smallest grains still have a 
finite geometrical cross section. 



